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INTRODUCTION 


This paper is a survey of computer-oriented power system analysis as practiced in the electric 
utility industry. Problems of the interconnected system in the western United States may be 
emphasized more than problems in other parts of the country because this is where the author 
gained most of his experience during eight years with the Bonneville Power Administration in 
Portland, Oregon. Of necessity, the survey can only highlight a few points. (An excellent source of 
background material is the July 1974 issue of the Proceedings of IEEE on computers in the power 
industry (ref. 1 ). 

Power systems are interconnected facilities for generating, transmitting, and distributing elec- 
tric energy. The backbone of a power system is a network of high-voltage overhead transmission 
lines that interconnect the powerplants with the load centers. Most lines are operated as balanced 
three-phase systems, but there are also high-voltage direct current links. Power systems usually 
extend over many states and comprise many utility companies. Figure 1 shows the major transmis- 
sion lines in the western United States and Canada, which are owned by about 40 utility companies 
who have formed the Western Systems Coordinating Council (WSCC) for joint studies on a com- 
bined systems basis. Joint studies are made by a technical staff located in Salt Lake City, which 
utilizes the computer facilities of the University of Utah. Similar “power pools” or “regional 
electric reliability councils” exist in other parts of the United States (fig. 2). 

In contrast to other forms of energy, there are few economic ways yet for storing electric 
energy on a large scale. One practical solution is the pumped-storage plant, where water is pumped 
into a reservoir at times of low load consumption and excess production at other plants, such as 
run-of-the-river hydroplants, and thermal plants running at minimum output. By and large, the 
energy demanded by the customers must therefore be generated at all times, following the load 
demand curve. Interconnections have been built primarily to take advantage of the greater diversity 
of load in a larger region to keep the total installed capacity as low as possible, and also to share 
reserve capacity in case of plant outages. Any imbalance between generation and utilization would 
require partial shedding of load through selective switching or voltage reduction. All generators on 
the system run in synchronism at 60 Hz in the United States and Canada, and any disturbance to 
the normal operation, such as an insulator flashover, leads to relative oscillations of the machines 
against one another. If these oscillations do not die out fast enough, then the system may become 
unstable and “collapse.” Therefore, the system must not only be balanced between generation and 
load under normal operation, but must also be designed to be stable against small and large 
disturbances. 

The growth of power systems is determined by the load demand, which has been doubling 
about every 1 0 years in the United States as well as in most other industrialized countries. Simple 
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arithmetic will show that such a geometric progression cannot continue forever. The Club of Rome 
and others deserve credit for drawing attention to the limits of growth. There is reason to believe, 
however, that the portion of electric energy in the total energy consumption, which is about 
25 percent today, may grow to about 50 percent by the year 2000, as more electric energy will be 
used for mass transportation, recycling, and other purposes. Therefore, electric power systems will 
continue to grow. 

A reliable supply of electric energy without imposed restrictions on consumption requires 
(a) large amounts of capital expenditure and (b) careful planning of the total system expansion 
through analysis of the normal operation as well as of the effects of disturbances. The word 
“expansion” must be emphasized because planning is almost always concered with additions and 
modifications of an already existing system. Power system planning is partly science and partly art 
and must answer questions of the type “where and at what time should a powerplant (transmission 
line) be built, what should be its rating, and should few units of large size be used or more units of 
smaller size?” For such studies, the load growth must be forecast in great detail for the next 5 years 
and in less detail for the next 20 years, always, of course, with a degree of uncertainty. Advance 
planning is essential because lead times (time from the decision to build something to its going into 
operation) are almost 5 to 10 years now for powerplants and 2 to 4 years for transmission lines. 

For those engaged in power system analysis and power system planning — the “software side” 
of power system research — it is well to remember that it is primarily the components of the system 
(such as powerplants, transmission lines, and circuit breakers) that require the most expensive 
research. Normally, it is not so much a matter of new technology but of continuity. 


HISTORIC PERSPECTIVE 


Systems analysis has always been important in the power industry simply because additions of 
powerplants and transmission lines require so much capital expenditure that their influence on the 
behavior of the overall system has to be analyzed before they are built. There is very little room for 
modifications after installation. There is reason to believe that power systems were the first systems 
studied in the modern sense of systems engineering. 

The mathematical foundations of power system analysis are very old. The well-known node 
and mesh equations were already explained in Maxwell’s books published in 1873 (ref. 2). He, in 
turn, relied on work done by Kirchhoff, Helmholtz, and others. It was recognized almost 50 years 
ago that meaningful system studies were impossible with hand calculations, which led to the 
development of specialized analog computers (“network analyzers”) in the 1920’s and 1930’s. 
Some are still in use today, and some are still being built for special purposes (e.g., for studies of 
electromagnetic transients). 


SOLUTION OF SPARSE NETWORK EQUATIONS 


Many problems in power system analysis lead to the solution of a system of linear equations. 
The steady-state behavior of an electric network is normally described by node equations of the 
form 
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m [ v] = m 


(i) 


where [T] is the nodal admittance matrix; [ V ] , the vector of voltages from node to datum; and 
[/] , currents injected into the node from datum. All matrices are normally complex. The diagonal 
element is the sum of all admittances of the branches connected to node i, and the off-diagonal 
element Y ^ is the negative admittance of the branch connecting nodes i and k. Since only few 
branches are connected to each node in power systems, [ Y] is “sparse” j(it has only a few nonzero 
entries). Nodes represent powerplants, load centers, and major substations, while branches represent 
lines, cables, and transformers. 

* < 

Gauss elimination 1 has clearly become the preferred method for solving linear equations' of 
the form of equation (1). It requires less memory and time than the Gauss-Jordan diagonalization 
process. In power system analysis, the elimination process for matrix [ K] is normally separated 
from the elimination process for the right-hand side [/] . This offers advantages if the system has to 
be solved repeatedly with the same matrix [Y] but with different right-hand sides [I], which is 
frequently the case. For such “repeat solutions,” only the process for the right-hand sides must be 
repeated. Often, [ Y] is symmetric; in that case, only the upper triangular matrix must be stored. 

All realistic power system problems are under the curse of dimensionality. Power flow studies 
must routinely by performed for networks of more than 1000 nodes with more than 2000 
branches. The major breakthrough in the solution of equations for such large power systems came 
in the early 1960’s with the exploitation of sparsity by ordered elimination (refs. 3,4). This has 
reduced solution times and memory requirements by almost a factor of 100 in smaller systems and 
much more in larger systems (see the ratio in fig. 3). It is believed that sparsity techniques were 
pioneered in the power industry, but mathematicians are well aware of it now (refs. 5—7). The 
following table illustrates the savings for a moderately sized system: 

Matrix data: 


Number of nodes 267 

Number of branches 423 

Number of nonzero elements above diagonal 

after elimination in 267 X 267 matrix 1015 

Solution times (IBM 7040): 

Triangular factorization of complex matrix 16.5sec 

Repeat solution 1.6 sec . 


Without sparsity, the upper triangular matrix would have approximately 36,000 elements, which is 
approximately 36 times more than in the table. Figure 3 (from ref. 8) compares the numerical 
effort of straight-forward matrix inversion with that of ordered elimination with exploitation of 
sparsity for typical power network problems. Figure 4 (also from ref. 8) shows the network graph 
and the sparse matrix after triangularization for a 62-node water distribution network, which 
illustrates that sparsity can also be exploited in non-electric problems. Another example would be 
the analysis of pin-jointed mechanical structures such as transmission towers (ref. 9). 

'Also called Gauss-Banachiewicz, triangulation, triangularization, triangular factorization, LU decomposi- 
tion, Gauss-Doolittle, Crout, Cholesky, etc. (sometimes in modified forms). 
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FAULT STUDIES 


Fault studies are made to Find the fault currents if faults occur at various locations in the 
system. Most faults involve only one phase of the three phases, for example, a flashover across an 
insulator to the tower. Fault current values are needed to check whether they lie within the 
interrupting capacity of those circuit breakers that remove the faulted line selectively and tem- 
porarily from the system. 2 They are also needed to set the sensing devices for the circuit breaker 
tripping mechanism . 

Classical fault studies solve the steady-state component of the fault current only, which is 
sufficient for the purposes mentioned before. The problem formulation leads to a system of linear 
equations of the form of equation (1). For balanced three-phase faults, a single-phase equivalent 
representation is used. For example, the voltage drop along the three phases of a line 



(where Z s is self impedance and Z m is mutual impedance) can be simplified for balanced conditions 
(lb = I a e'i l 20 °,I C = I a ei x 2 °°, analogous for V) to the single-phase equivalent 

dV a 

—fr-Vs-ZmVa- ( 3 ) 

The solution for phases B and C is the same as that for phase A, except for phasor rotations with a 
factor e~l 120 . The expression ( Z s ~ Z m ) is the equivalent impedance for balanced three-phase 
operation. 

For unbalanced faults, such as flashover from one conductor to the tower, symmetrical com- 
ponents are used. This is a well-established technique in power system analysis (ref. 1 0) whereby the 
three coupled phase equations such as equation (2) are transformed into three decoupled equations: 



2 The removal time (“dead time”) is kept as short as possible (just long enough to permit arc extinction in 
case of flashovers); it is noticed by the consumer only as a brief flicker. 
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with the linear transformation 
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for V and one that is identical for /, with a — e^ 120 ° . The new variables 0,1,2 are called zero, 
positive, and negative sequence components. With symmetrical components, one has to solve 
N equations three times rather than 3 N equations once, which saves memory and computer time. 
With sparsity techniques, however, the savings may no longer be as impressive because solution 
times increase about linearly with N rather than with N 3 . 


The techniques for classical fault studies are well developed, and it is unlikely that worthwhile 
improvements can be made. 


POWER FLOW STUDIES 


The flow of electric energy in interconnected ac systems is determined by the branch imped- 
ances, as expressed in Kirchhoff s laws, and cannot be controlled on an individual line. (There are 
exceptions, such as real power control with phase shifting transformers and high-voltage direct cur- 
rent links.) When lines are added to an existing system or when the effect of line outages is to be 
studied, it is therefore necessary to study the flows in the entire system. This is the power flow or 
load flow problem. Normally, only balanced conditions are studied with single-phase equivalents 
such as equation (3). Formulating the power flow equations in an iV-node system leads to AH nodal 
equations of the form of equation (1), except that the current is a function of the voltages, 


4 = 


p k iQk 
v k 


( 6 ) 


where is the real power into node k, Q ^ is the reactive power into node A:, and is the 
conjugate complex voltage from node k to ground. Normally, P ^ and Q ^ are specified, or P ^ and 
IK^I are specified. Equation (6) not only makes the power flow equations nonlinear but also 
nonanalytic because of the conjugate complex term. Therefore, the complex derivative is not 
defined, and numerical techniques based on derivatives must use pairs of real equations rather than 
complex equations and rectangular or polar coordinates rather than complex variables. 

The standard solution technique for power flow problems is now Newton’s method (ref. 1 1 ). 
The system of power flow equations 


[<?([*])]= 0 , 


(7) 
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with the unknown vector [x] , is solved iteratively with the system of linearized equations: 


and 



[Ax] =-[g] 


(8a) 


[*<*>] = [*<*-■>] + [Ax] . 


(8b) 


The Jacobian matrix [3g/3x] and the right-hand side -[g] in equation (8a) is evaluated at the 
approximate solution point [x^ ^] . Newton’s method only became practical for large power 
systems after sparsity techniques had been developed. The Jacobian matrix shows basically the same 
sparsity pattern as the admittance matrix in equation (1). Typically, a solution with good accuracy 
is reached in three to four iteration steps, independent of the size of the system. 


The Jacobian matrix gives a linearized model of the power flow equations around the solution 
point; therefore, it is very easy to calculate first-order sensitivities with repeat solutions, provided 
the triangularized Jacobian matrix has been stored. The equations are simply 

with [p] being those parameters for which the influence on the solution vector [x] is sought. 


The Jacobian matrix can also be used to calculate reduced gradients as the basis of optimiza- 
tion studies (refs. 12,13). (These techniques are beyond the scope of this paper.) 


STABILITY STUDIES 


Large disturbances, such as short circuits or powerplant outages, cause electromechanical 
transients in the form of relative oscillations between synchronous machines. These oscillations may 
be large enough to cause loss of synchronism in one or more machines. A generator that loses 
synchronism because of some disturbance is automatically disconnected from the system to avoid 
overheating and damage. Often, this increases the severity of the disturbance for other generators 
and, in turn, more generators may lose synchronism (“cascading outages”). This is the typical 
course of events in “blackouts.” 

Stability probems are more critical in geographically large systems (as in fig. 1 ) than in tightly 
meshed systems. Stability first became important in the 1930’s when hydroelectric plants were built 
far away from the load centers. At that time, “swing curves” (= rotor oscillations as a function of 
time relative to synchronous speed) were calculated to the crest of the first swing because power 
systems had enough damping that synchronism was practically never lost on subsequent swings. 
Today, stability must be checked beyond the first swing because modern fast-acting excitation 
systems on generators have decreased the system damping. The interconnection of formerly discon- 
nected power systems has also created new stability problems in the form of spontaneous oscilla- 
tions, which are sometimes damped with supplementary control signals in exciters and turbine 
governors. 
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Lyapunov’s second method has been proposed to find the region of stability directly without 
simulation. Many papers have been written on the subject because it is a challenging theoretical 
concept for scientific reasoning, but it is not yet a competitive alternative to the existing transient 
stability simulation programs. The region of stability obtained by Lyapunov’s method is too con- 
servative since the condition for stability is only sufficient but not necessary. In other words, the 
system might still be stable outside that region. The shortcomings of the method are overwhelming 
at this time and it seems questionable whether application to realistic power systems will ever be 
feasible (ref. 14), even though others disagree (ref. 15). 

Since there are no practical methods yet for assessing stability directly, simulating the behavior 
as a function of time, for specific disturbances assumed by the planner, remains the only practical 
alternative. Its main drawback is that the question of system stability is only answered for the 
specific disturbance, starting from specific initial conditions. Production-type stability programs can 
solve systems with up to 2000 nodes and 600 generators step by step with about 10 to 20 different 
exciter models and 5 to 1 5 different turbine governor models. 

In simulating electromechanical transients, two systems of equations must be solved simul- 
taneously, namely, the system of power flow equations. 3 

[£([*], Lv])] =0 (10) 


and a system of differential equations, 


dy 

dt 


[/([*],[y],0] 


(ID 


which describe the dynamic behavior of the turbine-generator rotors and of the exciters and turbine 
governors. 


Step-by-step solution methods for stability programs are classified in reference 16. Most pro- 
grams solve the differential equations and the power flow equations alternatingly, using a prediction 
of power flow state variables to solve the differential equations over one time step and using these 
results to obtain the new power flow solution at the end of the time step. There is a large spread in 
the magnitude of the eigenvalues in equation (11), assuming that the equations are linear or linear- 
ized of the form [ dy/dt ] = [A] [y] . Explicit methods such as fourth-order Runge-Kutta are very 
slow for such “stiff systems”; they require a small time step, dictated by the smallest time con- 
stants, which one would not anticipate from the smooth curves for the rotor oscillations. This 
“small time constant barrier” is being overcome with implicit integration schemes, such as the 
trapezoidal rule (ref. 16), which has been used quite successfully for electromagnetic transients in 
power systems since at least 1 96 1 (see ref. 1 1 in ref. 1 7). 


3 The rotor oscillations are slow enough to permit the use of steady-state equations for the network part. 


SUMMARY 


Only some topics of power system analysis could be described. While they cover the “big 
three” problems of fault, power flow, and stability studies, there are many other analysis problems 
in power systems which are partly summarized in reference 1 . 


One topic omitted here, the computation of electromagnetic transients, would have provided a 
good example of transfer of knowledge from one discipline to another, which is at least partly the 
objective of this workshop. The method of characteristics described in reference 17, which has 
become the standard solution technique for traveling wave problems on transmission lines, was first 
used to study pressure waves in hydraulic systems in the late 1920’s (see refs. 7 and 8 in ref. 17). It 
also illustrates that older techniques developed for hand calculations may still be valuable in our 
days of powerful computers, an observation which is also true for implicit integration with the 
trapezoidal rule. 

The solution techniques for most power system problems are highly developed by now and 
relatively efficient because of the exploitation of sparsity. Improvements are still possible, of 
course. 
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Figure 1. -Major transmission lines in the western United States. Source: Annual Report 1972 of 

the Western Systems Coordinating Council. 
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NATIONAL ELECTRIC RELIABILITY COUNCIL 



MARCA SUMMER 12,452 - WINTER 11,543 

(DUAL COORDINATION: MARCA & MAIN) 

MAIN SUMMER 27,085 - WINTER 23.546 

ECAR SUMMER 46,695 - WINTER 47,302 

NPCC SUMMER 39,606 
WINTER 45,212 


MAAC SUMMER 27,852 
WINTER 25,095 


SERC SUMMER 62,956 
WINTER 64,723 


SPP SUMMER 27,552 - WINTER 21,361 
(DUAL COORDINATION: SPP & MAIN) 

ERCOT SUMMER 20,252 - WINTER 13,693 


WESTERN SYSTEMS COORDINATING COUNCIL 


SUMMER 

WINTER 


56,470 (ACTUAL) 
60,180 (ACTUAL) 


Figure 2.— Regional electric reliability councils. Peak loads in MW for Summer 1972 and Winter 
1972/73 (forecast). Source: Sept. 1972 report of the National Electric Reliability Council. 



Figure 3.— Comparison of the numerical effort required by matrix inversion with that of ordered 
triangular factorization (OTF) for typical power network problems. 
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(a) NETWORK GRAPH. DRAWN (b) NONZERO PATTERN 

APPROXIMATELY TO O RESULTING MATRIX 

GEOGRAPHICAL A FILL-IN UPON TRIANGULARIZATION 

PROPORTIONS. 


Figure 4.— Sample small 62-node water distribution network. 
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